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^46 initio Molecular Dynamics Study of Glycine Intramolecular Proton Transfer in 

Water 

Kevin Leung and Susan B. Rcmpe 
Sandia National Laboratories, MS 1415 & 0310, Albuquerque, NM 87185 
(Dated: February 2, 2008) 

We use ab initio molecular dynamics simulations to quantify structural and thermodynamic prop- 
erties of a model proton transfer reaction that converts a neutral glycine molecule, stable in the gas 
phase, to the zwitterion that predominates in aqueous solution. We compute the potential of mean 
force associated with the direct intramolecular proton transfer event in glycine. Structural analyses 
show that the average hydration number (N w ) of glycine is not constant along the reaction coor- 
dinate, but rather progresses from N w =5 in the neutral molecule to N w =8 for the zwitterion. We 
report the free energy difference between the neutral and charged glycine molecules, and the free 
energy barrier to proton transfer. Finally, we identify approximations inherent in our method and 
estimate corresponding corrections to our reported thermodynamic predictions. 



I. INTRODUCTION 

Proton transfer plays diverse but fundamentally im- 
portant roles in aqueous solution phase chemistry. Trans- 
fer reactions in water, aided by the indistinguishable na- 
ture of protons, produce highly mobile charge carriers 
utilized in fuel cellsi and biological systems^ for electri- 
cal energy conduction. Protonation and deprotonation 
reactions on ionizable substrates activate and modulate 
biochemical reactions. For example, in Rubisco, the en- 
zyme that catalyzes the first reaction in carbon fixation, 
many steps in the reaction mechanism consist of proto- 
nation or deprotonation events on the enzyme-substrate 
complex. 3 In biological ion channels, changes in protona- 
tion states along the ion conduction path modify trans- 
port function^ and changes in protonation state in the 
gating region of some channels presumably trigger open- 
ing or closing of the channels^ In order to understand 
such chemical processes associated with proton transfer 
in molecular detail, it is crucial to have methods that 
yield accurate estimates of the accompanying free energy 
changes. 

Free energies, which characterize molecular interac- 
tions and govern the likelihood and rates of chemical re- 
actions and conformational changes, are among the most 
important properties to be computed in aqueous envi- 
ronments. Since its introduction, the ab initio molecular 
dynamics approach (AIMD) has demonstrated remark- 
able success in predicting and modeling the hydration 
structure, vibrational signature, and electronic proper- 
ties of ions and molecules in water & In recent years, it 
has also found increasing utility in computing thermo- 
dynamic properties^ including free energy differences by 
potential of mean force methods. 8 ^ These methods have 
successfully reproduced several pKa's for deprotonation 
reactions in water. 

In fact, the recent introduction of various methods 
to calculate free energies, such as reaction-path finding 
techniques fiS targeted molecular dynamics^ and related 



methods, and a new theory of the liquid state that em- 
phasizes the local hydration environment^ has created 
the potential for combining AIMD and unbiased, rigor- 
ous predictions of free energy barriers and differences. 
Relatively few aqueous phase free energy calculations ex- 
ist in the literature based explicitly on AIMD data. To 
benchmark new advances in free energy methods, AIMD 
free energy calculations on simple, representative systems 
(especially proton transfer reactions) would be extremely 
valuable. 

A convenient and popular model system for studying 
proton transfer events in aqueous environments is the in- 
tramolecular proton transfer of glycine in water. Water 
triggers zwitterion formation from the neutral molecule 
stable in the gas phase. Since the stability of the zwitte- 
rion relative to neutral conformers in water is entirely due 
to glycine- water interactions, 1 " theoretical predictions for 
glycine in water are sensitive to the parameterization of 
water-glycine functional group interactions. This makes 
the glycine intramolecular proton transfer reaction a sen- 
sitive and valuable benchmark for calculating hydration 
effects on biological functional groups. In fact, when new 
methods are devised to sample molecular conformations 
at finite temperature quantum mechanically, glycine is 
often a model of choice as a first application!^^ 

Part of the reason for the popularity of glycine as 
a model for proton transfer studies in proteins is that 
its tautomerization free energy difference and barriers 
have been measuredii£iiLi2i±2i22i2ii22 Glycine can exist 
in many neutral form conformers in addition to the zwit- 
terion tautomer, as illustrated in Fig. 1. In the gas phase, 
the zwitterion (ZW) is unstable and spontaneously col- 
lapses to the neutral form (NF) conformer "Hp. "^24 In 
contrast, the ZW is more stable than any neutral form 
in water due to large local charges and a large dipole 
moment. A free energy barrier of 14.3 kcal/mol is asso- 
ciated with the ZW^NF interconversion^SS and ZW is 
more stable by 7.27 kcal/molii£ From these results, the 
reverse, NF^ZW reaction is estimated to have a free en- 
ergy barrier of ~ 7 kcal/mol and there are indications 
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FIG. 1: Glycine tautomers. Black circles: H; grey: C; small 
and large open circles: N and O, respectively. Panel (a): zwit- 
terion (ZW) glycine; (b)-(d) three stable/metastable neutral 
form (NF) glycine molecules with configurations optimized in 
the gas phase, respectively, Hp, Ip, and IIIpi— Note that some 
works in the literature label 0(1 ) and 0(2) in the reverse way. 

that the ZW— >NF free energy barrier is mainly entropic 
in origin— i Other experimental data provides additional 
information on glycine hydration. Neutron scattering re- 
sults on the hydration structure of the ammonium group 
in concentrated glycine solutions have been reportedpS, 
as well as mass and size-selected photoelectron spectro- 
scopic studies of hydrated glycine anions^ 

While experimental data makes glycine a choice sys- 
tem for study, comparison of theoretical and experimen- 
tal studies of the proton transfer between the ZW and 
NF tautomers is complicated by the fact that the most 
stable aqueous NF conformer has not been identified^ 
The various NF conformers are related by rotation of the 
C - C, C - O, and C - N bonds. Conformer Ip (Fig. 1(c)) 
is the most stable structure in the gas phase while Hp is 
~ 1 kcal/mol higher in energyi£L2£ Intuitively, a direct 
intramolecular proton transfer from ZW to NF in wa- 
ter should go first through Hp, which is closest in struc- 
ture to ZW among the neutral form conformers. This 
is the proton transfer reaction that is most accessible 
to finite temperature, quantum mechanic (or quantum 
mechanics-based) simulations. Thus practitioners using a 
molecular dynamics simulation, based on a reactive force 
field fitted to Hartree-Fock data^?' 3 - have claimed to re- 
produce the experimental values for ZW— >NF free en- 
ergy difference and barrier by just considering the Hp 
neutral form. Other quantum chemistry works, how- 
ever, find a small reverse IIp^ZW barrier and interpret 
Hp as an intermediate, not the final or most stable NF 
product>iiii22*2i To further complicate interpretation of 



the results, these quantum chemistry results depend on 
the details of the calculations and so are not necessarily 
in quantitative agreement with each other. 

In this work, we do not seek the most favorable form 
of (metastable) NF glycine in aqueous solution. Instead, 
we use the ab initio molecular dynamics (AIMD) ap- 
proach to focus on the direct IIp^ZW proton trans- 
fer in water. AIMD treats the valence electrons of all 
atoms quantum mechanically, using density functional 
theory (DFT). It also samples glycine and water confor- 
mations at finite temperature via its molecular dynam- 
ics capability. While AIMD is computationally costly, 
it has several important advantages when treating in- 
tramolecular proton transfer in glycine, (a) Glycine it- 
self is treated quantum mechanically, which is pertinent 
to the breaking and making of chemical (covalent) bonds 
that take place during proton transfer, (b) The num- 
ber and conformations of water molecules in the first hy- 
dration shell, and the hydrogen bonds they form with 
the glycine atoms, are allowed to fluctuate and vary. 
This is significant because hydrogen bond network fluc- 
tuations have been known to be crucial in other proto- 
typical proton transfer reactions in water i22 (c) These 
first hydration shell water molecules are treated quan- 
tum mechanically. Widely used force field models for 
water and the carboxylate (-COO - ) functional group, 
which is a crucial part of ZW glycine, predict carboxylate 
group hydration structures^i^iS^SiSLS that disagree 
with experiments^24£ This is true even if the carboxylate 
group containing ion is treated quantum mechanically 
and only water is treated classically— (i.e., a QM/MM 
treatment). In contrast, AIMD predictions of hydra- 
tion structure are in good agreement with experiments. 43 
Likewise, in the one case where a QM/MM estimate of 
the ZW/IIp free energy difference is reported, the the- 
oretical prediction overestimates the experimental value 
by a factor of four 44 This emphasizes the sensitivity of 
this free energy difference to glycine-water interaction pa- 
rameters. 

The role of water conformational changes and water- 
glycine hydrogen bond fluctuations on the glycine in- 
tramolecular proton transfer reaction is central to this 
work. We dynamically sample water conformations us- 
ing AIMD. These water configurational changes may be 
especially important for the ZW— >IIp reaction. This is 
because the glycine dipole moment decreases along the 
reaction coordinate, leading to a reduction in the mean 
hydration number. Our work will reveal trends concern- 
ing hydration number variations and hydrogen bond net- 
work fluctuations as the intramolecular proton transfer 
proceeds. 

A number of pioneering quantum chemistry works on 
glycine tautomerization have used a quantum treatment 
of the glycine molecule as well as a few water molecules in 
the glycine first hydration shell (i.e., a "supermolecule" 
approach). A static optimization of the supermolecule 
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geometry is performed, and a polarizable dielectric con- 
tinuum is used to treat the outlying waterji^i^ The ad- 
vantage of such methods is their relative computational 
ease, which allows the sampling of many glycine neutral 
form conformations. 

In connection with these studies, our AIMD predic- 
tion of the progression in the average hydration num- 
ber around each glycine atom, as the reaction coordinate 
varies, may contribute to future supermolecule studies of 
proton transfer reactions in water. Working within the 
framework of a statistical mechanical theory of liquids 
called quasi-chemical theory^ for example, accurate hy- 
dration free energies for metal ions have been computed 
by forming supermolecular clusters containing the metal 
ion and the full set of water molecules in the first hydra- 
tion shell and embedded in a dielectric continuum model 
of bulk water Clusters containing one glycine 
and 3 or 6 water molecules have been considered in re- 
cent published quantum chemistry studiespi^i^ As will 
be shown in this work, this is not sufficient to complete 
the first hydration shell of ZW glycine. We will make a 
preliminary attempt to quantify this effect. 

The paper is organized as follows. Section 2 describes 
the method used. Our predicted potential of mean force 
and correlation functions are described in Sec. 3. Sec- 
tion 4 concludes the paper with further discussions. In 
the appendices, we discuss corrections to the AIMD re- 
sults by estimating the effect of zero point energy contri- 
butions and of using different exchange-correlation func- 
tionals, and also we address the effect of using a finite- 
sized simulation cell. 



II. METHOD 

We perform ab initio molecular dynamics simulations 
on a system with 1 glycine and 52 H2O molecules. Finite 
size effects will be addressed in Appendix B using a sim- 
ulation cell with 98 H2O molecules. The Car-Parrincllo 
Molecular Dynamics (CPMD)£i code is applied, along 
with the BLYP gradient corrected exchange correlation 
functional£2i£i and Troullicr-Martins pseudopotcntials. 54 
BLYP has been shown to yield water-water pair cor- 
relation functions g(r)i as well as hydration structures 
of NH^ and HCOO~ ions that are in good agreement 
with experiments^ 3 *^ The effect of using other, perhaps 
more accurate, exchange correlation functionals will be 
addressed in Appendix A. 

The simulation box is cubic with linear dimension 
11.76 A, which corresponds to a water density of 1.00 
g/cm 3 plus the experimental glycine zwitterion volume 
of 72 A 3 ^ The time step used is 5 a.u. (0.121 fs), and 
the deuterium mass is assumed for all protons through- 
out, although we will continue to use the word "proton." 
The temperature is kept at T=300 K using a thermostat 
unless otherwise stated. 



We compute the potential of mean force, AG(R), 57 
along a reaction coordinate defined as the difference be- 
tween the N - H and 0(1) - H distance: 

R = i?N-H - -Ro(l)-H- (1) 

0(1) is the Hp glycine oxygen atom covalently bonded 
to the acidic proton, and H is the proton being directly 
transferred between 0(1) and the nitrogen atom. Pro- 
gression along this coordinate represents a direct proton 
transfer from the ZW to the Hp form of glycine. For 
the Hp neutral conformer (Fig. lb), i?o(i)-H ~ 1 A, 
and N and H are not covalently bonded, making -Rn-h 
larger than 1.4 A, and hence R > 0. For the ZW glycine, 
-Rn-h ~ 1 A, 0(1) - H is not covalently bonded, and 
R < 0. We use the umbrella sampling method^ (with 
harmonic constraints or biasing potentials) to compute 
the free energy profile along this reaction coordinate. 
Since one of the sampling windows is the unconstrained 
glycine zwitterion, this allows us to examine the structure 
of water around stable, equilibrium zwitterion glycine, 
and compute the pertinent water-glycine pair correlation 
functions (g(r)) as welli^ 

The barrier height AG* depends on the reaction path- 
way. A more sophisticated approach would be to use the 
path-sampling method^ which in principle can sample 
the possible reaction pathways weighted by the proper 
statistical mechanical weight. This is computationally 
expensive, however, and seldom implemented within an 
AIMD setting. For the relatively simple system of glycine 
intramolecular proton transfer, we expect our relatively 
simple reaction coordinate to be adequate. This is be- 
cause the reaction coordinate is similar to the one con- 
sidered in Ref. 31 - it describes a straightforward direct 
proton transfer event. Using a simple N - H bond dis- 
tance reaction coordinate, Kassab et al. obtain results 
similar to their own gas-phase transition state finding 
calculations. 31 Furthermore, when we initiate a glycine 
Hp conformation in water, we observe that it undergoes 
a direct, intramolecular proton transfer to the ZW form 
in less than 1 ps. This observation corroborates the re- 
action coordinate we have chosen. 

While the two oxygen atoms in ZW are in principle 
equivalent, they interconvert on a relatively long time 
scale, and they can be treated as distinct for the dura- 
tion of our AIMD trajectories. In contrast, we do ob- 
serve rotation around the N - C bond in our simulations, 
which is also known to occur in picosecond timescales. 61 
The three protons bonded to the nitrogen atom in the 
ZW glycine are therefore equivalent in our simulations, 
and, at any time step, the relevant proton in Eq. 
is taken as the one closest to 0(1). 59 To compute the 
free energy profile along R, a harmonic penalty function 
V(R) = Ai(R~Ri) 2 is used to constrain R in six umbrella 
sampling windows with a progression of target distances 
where Ai ranges from 0.7 eV to 1.2 eV/AA (see Ta- 
ble^for details). Window (a) corresponds to ZW glycine, 
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window 


a (ZW) 


b 


c 


d 


e 


f (Hp) 


Ri (A) 


0.0 


-0.8 


-0.5 


-0.1 


0.2 


0.7 


Ai (eV) 


0.0 


1.0 


1.0 


1.2 


1.0 


0.7 



TABLE I: Constraint parameters A and B for the umbrella 
sampling windows. See Eq. 



which is stable and requires no constraint. 

AIMD trajectories are initiated by inserting a ZW or 
Hp glycine molecule into a simulation box with 54 wa- 
ter molecules, which has been equilibrated previously us- 
ing the empirical SPC force field for water Two water 
molecules overlapping the glycine molecule are removed. 
We first conduct a 10 ps QM/MM molecular dynamics 
trajectory for ZW and Hp glycine.™ The resulting molec- 
ular configurations are used as the starting points for the 
AIMD trajectory in the ZW and Hp windows (i.e., win- 
dows (a) and (f)). Using AIMD, the QM/MM Hp con- 
figuration is equilibrated for 4 ps, and then statistics are 
collected for the next 10 ps. The window (f) constraint 
(Table P) is always applied to keep the Hp glycine from 
undergoing intramolecular proton transfer. The starting 
configuration of window (e) is taken 8 ps into the Hp tra- 
jectory, equilibrated at the new constraint (Table for 
6 ps, and run for another 10 ps. Window (d) is spawned 
from window (e) 5 ps into its trajectory, then equilibrated 
for 2 ps with the new constraint. As will be discussed, in 
windows (d) and (e), statistics are collected for 20 ps. 

As for window (a), we first use AIMD to re-equilibrate 
a QM/MM ZW configuration at T=300 K for 10 ps. 
The final configuration of this trajectory is once again 
re-equilibrated at T=350 K for 1 ps, and then statistics 
are collected at this temperature for 10 ps. Both win- 
dows (b) and (c) are spawned from the configuration at 
the end of the window (a) trajectory, equilibrated for 
1 ps with the new constraint, and then statistics are col- 
lected for 10 ps. These two windows are thermostated at 
T=350 and 300 K, respectively. The reason for using a 
higher temperature, and the small effect this has on the 
potential of mean force, will be discussed later on. 

Our AIMD trajectories within each sampling window 
last 10 or 20 ps. Our total AIMD trajectory length is 
similar to that used to correctly reproduce experimen- 
tal results for deprotonation of histidine residues. 8 Note 
that the reactive force field based glycine intramolecular 
proton transfer work of Ref. 28 uses trajectories only a 
few times longer than the one reported in our work. Fur- 
thermore, our use of longer trajectories in the windows 
that contribute most to AG, and a higher temperature 
for the ZW window that only weakly contributes to AG 
(see below) further reduces statistical uncertainties. We 
estimate the cumulative statistical uncertainty of our pre- 
dicted AG and AG* between Hp and ZW glycine to be 
of order ~ 1.4 kcal/mol. (AG* for the reverse reaction, 
Hp— >ZW, is much smaller.) As a result, one limitation 
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FIG. 2: Sample hydration structures around just the car- 
boxylate/carboxylate acid group (not the amine group) in 
aqueous glycine. Only the first hydration shell is shown. Th 
grey, blue, red, and white spheres represent C, N, O, and H 
atoms, respectively. Upper left panel: long lived zwitterion 
hydration shell configuration at T=300 K. Upper right: zwit- 
terion hydration shell after a trajectory of 5 ps at T=350 K. 
Lower left: window (d) (transition state region); lower right: 
window (f) (neutral form glycine). 



of our AIMD work is that we cannot resolve the temper- 
ature dependence of AG for this intramolecular proton 
transfer reaction. Nevertheless, the potential of mean 
force computed using AIMD potentially can be used to 
calibrate or refine force fields, which then can be applied 
to address the entropic contribution to AG accurately 
and efficiently. 



III. AB INITIO MOLECULAR DYNAMICS 
RESULTS 

In this section, we first consider the hydration struc- 
tures of the glycine functional groups. We also inves- 
tigate the time dependence of the hydration numbers to 
demonstrate that the simulation time we use in each um- 
brella sampling window is adequate. Then we report the 
potential of mean force and hydration structures along 
the ZW— dip reaction coordinate. Finally, we discuss our 
results in relation with experiments and quantum chem- 
istry based supermolecule calculations. 



A. Zwitterion hydration structure 

ZW glycine, the most stable tautomer in water, has 
a large gas phase dipole moment and charge separation 
along the molecular framework that leads to strong in- 
teractions with water. Despite this, we expect, on the 
basis of work with formate ion hydration^ that the co- 
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FIG. 3: Time dependent hydration numbers for O(l) and 
0(2) combined, (a) Window (a); (b) window (c); (c) window 
(d). Based on the first minimum of the glycine- water g(r), we 
define a water molecule to be in the hydration shell if it one of 
its protons is within 2.5 A of the glycine oxygens. Filled and 
hollow triangles indicate where water molecules enter the hy- 
dration shell of O(l) or 0(2) for the first time, or leave those 
hydration shells for at least 1 ps. There are numerous tran- 
sient fluctuations and re-entrances into hydration shells. For 
concreteness, we arbitarily define a water molecule as "first" 
entering the hydration if it resides there for at least 1 ps later 
in the trajectory, and if it was not previously in the hydration 
shell for at least 2 ps. 



ordination number of the carboxylate oxygens will expe- 
rience large fluctuations^ We indeed find that the com- 
bined instantaneous hydration number for these oxygens 
ranges from 3 to 7, as illustrated in Fig. 3. At T=300 K, 
however, the first hydration shell of the -COO - group ex- 
hibits relatively slow dynamics during the first ~ 10 ps of 
the AIMD trajectory. For instance, in a 5 ps stretch, one 
of the water molecules bound to 0(2) briefly leaves the 
hydration shell, is replaced by a second water, and then 
returns, displacing the newly added water molecule. This 
long-lived ZW configuration has two water molecules hy- 
drogen bonded to 0(1) and three to 0(2) (see Fig. 2). 

We did not observe similar slow dynamics in AIMD 
simulations of the aqueous formate ioni 3 . using the BLYP 
exchange correlation functional, despite the fact that 
HCOO - is a charged species. We speculate that, whereas 
the carbon end of the formate ion is hydrophobic, the 
-NH 3+ group in the ZW glycine forms another strong 
hydration shell and concentrates water molecules in the 
vicinity of the -COO - group, causing a more persis- 
tent carboxylate-water interaction in the ZW glycine. 



Regardless of the reason, while this long-lived hydra- 
tion structure persists, the reaction coordinate R will be 
locked into a relatively narrow distribution. As a result, 
sufficient statistics for the ZW glycine cannot be accu- 
mulated using a 10-20 ps trajectory. Instead, we collect 
statistics for windows (a) and (b) at T=350 K. The dy- 
namics at this higher temperature are much faster and 
involve several exchanges of water molecules between the 
first hydration shell and the bulk liquid. See Fig. 3. 

As R increases, the glycine dipole moment decreases, 
its interaction with water weakens, and the dynamics 
of water molecules around the glycine molecule become 
faster. Figure 3(b) shows that the dynamics in window 
(c) and (d) are comparable to that in (a), despite the fact 
they are run at T=300 K. This indicates that an elevated 
temperature is not needed for window (c) for a 10 ps tra- 
jectory. Figure 3 actually depicts two complementary 
quantities, the instantaneous hydration number of 0(1) 
plus 0(2), and the times at which the water molecules 
first enter or leave the hydration shells. The former in- 
cludes all transient fluctuations and the rapid motion due 
to water molecules briefly forming and breaking hydrogen 
bonds with 0(1) or 0(2). The latter filters out the rapid 
fluctuations, but its definition is somewhat arbitary; we 
require that water molecules have at least a 1 ps residence 
time to be so counted. At even higher temperatures, such 
that the time scale of exchange between first hydration 
shells and the bulk water region is faster than 1 ps, this 
criterion will have to be redefined. 

We note that the residence times of water molecules in 
the hydration shells are sensitive to the strength of inter- 
action between carboxylate oxygens and water molecules; 
a slight increase of just 1 kcal/mol will increase them 
significantly. Hence such quantities are expected to be 
sensitive to temperature as well as to the exchange cor- 
relation functional used4£ 

Figure 4 depicts the ZW glycine-water pair correla- 
tion functions, g(r), for selected atomic sites. The two 
carboxylate oxygen atoms form an average of 4.7 total 
hydrogen bonds with water molecules. 0(2) exhibits a 
hydration number N w ~ 2.66, in good agreement with 
the hydration number of the formate ion computed using 
AIMD4 3 - 0(1) forms an average of 2.0 hydrogen bonds 
with water molecules. It also forms an intramolecular 
hydrogen bond with one of the ammonium group pro- 
tons 40% of the time, assuming a hydration shell ra- 
dius of 2.5 A. The combined hydration number of 4.7 
is smaller than that of 5.3 predicted by Alagona et alm^ 
using empirical force fields, and is much smaller than the 
hydration number of 7 for carboxylate groups predicted 
by Jorgensen and Gao's OPLS force fields. 38 Integrating 
g(r) between the ammonium protons and water oxygen 
atoms up to its first minimum, the ammonium group is 
found to have 3.0 water molecules in its first hydration 
shell. Using instead the g(r) between the nitrogen and 
water oxygen, we find N w » 4 under its first peak, with 
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FIG. 4: Pair correlation functions between four atoms in ZW 
glycine, and oxygen (dashed line) and hydrogen atoms (solid 
line) on water molecules, (a) O(l): (b) 0(2): (c) N; (d) one 
of the three equivalent protons on the ammonium group. N w 
indicates the number of water molecules around the atom, 
found by integrating to the first minimum in g(r) between 
the solute and water oxygen atoms. 

an average N - O distance of 2.75 A. A neutron scatter- 
ing study for 5 mol % glycine solution deduces a similar 
mean N - O distance of 2.85 ± 0.05 and N w = 3.0 ± 0.6, 
based on analysis of the N - O correlation^ Since the 
ammonium protons are directly hydrogen bonded to the 
water oxygens while the nitrogen atom is not, we deem 
the former more useful in determining the total number 
of hydrogen bonds, particularly because of the elevated 
temperature used in our simulations. Hence we estimate 
that AIMD predicts 3.0 hydrogen bonds between the am- 
monium group of ZW glycine and water molecules. 

Unlike the carboxylate oxygens, the combined coor- 
dination number of the three protons in the -NH3" do 
not exhibit large fluctuations, and thus these fluctuations 
should not contribute significantly to AG(R). 



B. Potential of mean force and hydration structure 

Figure 5 depicts AG(R) along the reaction coordinate 
R. The ZW glycine is found to be more stable than the 
Hp conformer by a free energy difference, AG, of 11.2 
kcal/mol. We also find a 12.7 kcal/mol transition state 
barrier, AG*, between these two forms of glycine located 
at R w 0.2 A. The statistical uncertainties in both AG 
and AG* are estimated to be about 1.4 kcal/mol. On the 
average, this transition state structure exhibits 0(1) - H 
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FIG. 5: Potential of mean force as a function of the cho- 
sen reaction coordinate R = Rn-h — Ro-h- The six um- 
brella sampling windows are demarcated with dashed lines, 
and snapshots of the glycine molecule in some of the windows 
are depicted. 

and N - H distances of ~ 1.1 and ~ 1.3 A, respectively. 
We will consider several corrections to AG and AG* be- 
low, but they will not qualitatively modify these conclu- 
sions. 

Figure 6 plots as functions of R the mean hydration 
numbers for 0(1) plus 0(2), N, and the proton be- 
ing transfered. They are compiled across all six um- 
brella sampling windows. The results at the boundary 
of two windows are averaged. As R increases, the ZW 
glycine continuously transforms to the Hp neutral form 
and the glycine dipole moment decreases, water becomes 
less structured around it, and the hydration number de- 
creases for all glycine atoms we examined. In particu- 
lar, when R reflects a Hp neutral form configuration, the 
NH2 and COOH groups do not strongly bind to water 
molecules; each amine proton forms hydrogen bonds to 
only an average of ~ 0.5 water molecules. As a result, 
the hydration numbers drop almost by a factor of two in 
these cases. 

As for the proton being transfered, it forms covalent 
and intramolecular hydrogen bonds with N and 0(1) in 
the range -0.8 A< R < 0.7 A, but forms no hydrogen 
bonds with water molecules. If R is further increased, 
the intramolecular hydrogen bond cannot be sustained, 
and the proton begins to establish an intermolecular hy- 
drogen bond with water molecules. Thus hydrogen bond 
network fluctuations are crucial to proton transfer, as ob- 
served in proton transport in water ml The statistical un- 
certainties in these i?-dependent hydration numbers are 
of order ~ 0.3. Considering this, the hydration numbers 
are reasonably well converged, smooth functions of R. 
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FIG. 6: independent mean number of water molecules in 
the first hydration shell, N w , for various atoms, (a) O(l) 
and 0(2) combined; (b) N; (c) H involved in intramolecular 
proton transfer. The cut-off criteria for determining these 
coordination numbers are: O - H w distance of 2.5 A; N - O w 
distance of 3.5 A; and H - O w distance of 2.4 A; respectively, 
where the subscript "w" refers to water molecules. 



The predicted progression of hydration numbers should 
be valuable for future quasi-chemical calculations. 46 

At the transition state, where R ~ 0.2 A, the hydration 
numbers (Fig. 6) closely resemble those found for the Hp 
neutral form conformer, but differ significantly from that 
of the ZW. 

To estimate the statistical uncertainties in our ther- 
modynamic quantities, we define A AG = |AG(i?2) — 
AG(R\)\. For windows (b)-(e), i?2 (-Ri) is the maximum 
(minimum) R used in each window in Fig. 5. For win- 
dows (a) and (f), Ri and R\ are the R values at the 
local minima, respectively. The statistical uncertainty 
in AAG is estimated by dividing the 10 ps trajectory 
in each window into 4 blocks of 2.5 ps each. For win- 
dows (d) and (e), we run a trajectory twice as long. The 
estimated deviation from the mean computed using 2.5 
and 5 ps blocks are almost identical for window (d), sug- 
gesting that the correlation time for AAG is less than 
2.5 ps. Recall that Fig. 3(c) shows only 5 instances each 
of "new" water molecules entering, and completely leav- 
ing, a hydration shell within the 10 ps trajectory segment 
depicted in this window. There are also numerous tran- 
sient fluctuations in N w , however, and these fluctuations 
allow statistically meaningful sampling of AG(R). The 
correlation time appears slightly larger than 2.5 ps for 
window (e). Given that most windows are sampled for 
10 ps, we report the cumulative uncertainty in AG as 



twice our estimated overall deviation from the mean es- 
timated using 2.5 ps time blocks. 

Recall also that windows (a) and (b) are sampled at 
T=350 K. For these two windows, we estimate AG(i?) 
at T=300 K by assuming that the free energy bar- 
rier is mainly due to entropyiSi Then AG(R)/ (k-gT) oc 
log P(R), where P(R) is the probability, assumed to be 
temperature independent, that the reaction coordinate 
spontaneously exhibits value R. As shown in Fig. 5, 
AAG in these two windows combined is 1.4 kcal/mol. 
Even if AG(R) and P(R) are entirely due to enthalpy, 
the resulting small error arising from sampling these win- 
dows at T=350 K will be a small fraction of 1.4 kcal/mol, 
and will have little effect on the ZW— dip AG. 



C. Neutral form glycines: spontaneous direct 
proton transfer 

A typical Hp neutral form glycine hydration structure 
is depicted in Fig. 2. The carboxyl oxygen and the two 
amine group protons of neutral form glycine exhibit hy- 
dration numbers of 2.1 and 1.0, respectively. These num- 
bers are obtained using the constraint of Table [I] 

In fact, we find that a constraint is necessary to stablizc 
the Hp glycine molecule in AIMD and QM/MM trajec- 
tories. Upon releasing the constraint, the Hp molecule 
spontaneously undergoes a direct proton transfer to the 
ZW form in a sub-picosecond time scale, regardless of 
the starting configuration. This is consistent with the 
small proton transfer free energy barrier computed using 
umbrella sampling. 



D. Comparison with experiments 

We predict a ZW— dip activation barrier (AG*) of 12.7 
kcal/mol, which is ostensibly within 1.6 kcal/mol of the 
reported experimental AG* ^ We note that the experi- 
mental rate has been measured with both nuclear mag- 
netic resonance 16 and the thermally modulated chemi- 
cal relaxation method^ and the reported intramolecular 
proton transfer rates are within 12% of each other. We 
also find a free energy difference for conversion from the 
ZW to the Hp conformer (AG) of 11.2 kcal/mol, which 
is ostensibly 54% larger than the experimental value of 
7.27 kcal/moliL 6 ^ As will be shown, various corrections 
to our AIMD free energies are small. 

Thus we find, in qualitative agreement with quantum 
chemistry calculations that the experimental results 
do not correspond to the neutral form Hp conformer that 
is the focus of this study. The NMR result states that 
the rate it measures is associated with a NF product 
that must undergo proton exchange with water before 
reverting back to the ZW form; the Hp conformer, which 
we predict to have a picosecond lifetime, will not qualify. 
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It is likely that the observed rateiS^ involves the ZW 
first tranforming into the Hp molecule, and then on to 
one or several more stable NF conformcrs that undergo 
proton exchange in water and lead to the coalescence of 
NMR line shapes^ The NMR free energy barrier will 
reflect a composite of these processes. 

We have not computed the free energies and barriers of 
other NF conformcrs compared to the Hp neutral form. 
From our discussion in the previous section, our predicted 
AG = 11.2 kcal/mol exceeds the experimental value of 
~ 7 kcal/mol, which is consistent with the fact that Hp is 
only an intermediate, and a lower free energy conformer 
exists. The search for this stable conformer using AIMD 
will be left to future work. 

We find that AG predicted by AIMD and the super- 
molecule approach, which uses one glycine and three wa- 
ter molecules plus dielectric continuum, differ by several 
kcal/mol if BLYP is used in the latter case. This com- 
parison will be presented in Appendix A along with the 
suggestion that the discrepancy exists because the super- 
molecule lacks a fully occupied first hydration shell. 



IV. CONCLUSIONS 

Using ab initio molecular dynamics calculations and 
the BLYP exchange correlation functional, we find a 
12.7 kcal/mol free energy barrier between zwitterion 
and conformer Up of the neutral form glycine in water. 
The statistical uncertainty is estimated to be of order 
1.4 kcal/mol. We predict a 11. 2± 1.4 kcal/mol free energy 
difference between the zwitterion and this Up conformer. 
The experimental free energy difference between the zwit- 
terion and the neutral form is 7.27 kcal/mol, although 
precisely which neutral form conformer dominates in wa- 
ter has not yet been determined experimentally, and will 
be addressed in future work. 

We also gain useful qualitative insight on hydration 
structures from our AIMD simulations. We find that the 
hydration structure of the -COO - group in the zwitterion 
is similar to that of the formate ion^& forming an average 
of 4.7 hydrogen bonds with water, while the -NHg" group 
forms 3.0 hydrogen bonds with water molecules, yielding 
8 water molecules in the nearest hydration shell of the 
ZW. The carboxyl oxygen and the two amine group pro- 
tons of neutral form glycine exhibit hydration numbers 
of 2.1 and 1.0, respectively. The coordination numbers 
along the reaction coordinate interpolate between these 
limits. 

This work demonstrates the viability of AIMD to pre- 
dict free energy changes in aqueous reactions. AIMD 
allows fairly extensive sampling of the first hydration 
shell water configurations, which is found to have strong 
correlation with the progress along the reaction coor- 
dinate. While our conclusions about the crucial role 
of neutral form glycine conformers other than Up is in 



agreement with some quantum chemistry supermolecule 
calculations (and in disagreement with reactive force 
field works 28,29 ), our quantitative results may be particu- 
larly valuable toward accurate parameterization of future 
quasi-chemistry calculations of proton transfer in aque- 
ous environments^ 
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Appendix A: Corrections to ab initio molecular 
dynamics results 

Several corrections to the thermodynamic quantities 
(AG and AG*) determined from our AIMD trajectories 
should be considered. We show that they will not alter 
our conclusions. These include zero point energy correc- 
tions and use of a more sophisticated exchange correla- 
tion functional. 

To estimate such corrections, we conduct gas phase 
cluster energy minimization calculations using the Gaus- 
sian code, with B3LYP — and BLYP exchange correla- 
tion functionals and various basis setsS Cluster opti- 
mizations are carried out in the 6-31G(d) basis, to re- 
produce results found in the literature, as well as in the 
6-31+G(d,p) basis. Frequency calculations confirm that 
the predicted structures are true minima, and zero point 
energies are computed at the same level of theory. Re- 
fined single point estimates of energies are obtained using 
an extended 6-311+G(2d,p) basis set applied to the con- 
figurations found in the minimization calculations. 

A. Zero Point Energy 

The results in Sec. lIIIBl do not include zero point en- 
ergy (ZPE) corrections. It is fairly costly to include ZPE 
effects in AIMD simulations via the path integral formal- 
ism. ZPE corrections, however, can be estimated during 
post processing of the trajectory data by appealing to 
supermolecule calculations. Here we consider a super- 
molecule geometry similar to that in Ref. 31. This cluster 
has one glycine and three water molecules, as shown in 
Fig. 7, and a polarizable dielectric continuum modelS is 
used to treat the bulk water boundary conditions. Using 
the basis and refined basis sets described above, we find 
that ZPE indeed raises the free energy of ZW glycine by 
only ~ 1 kcal/mol relative to Hp, while the transition 



9 







mS- 




(c) 


<d) © 












> 



FIG. 7: Structures of a gas phase cluster with one glycine plus 
3 water molecules, similar to Ref. 31. (a) Hp; (b) transition 
state; (c) ZW; (d) bridging water. 



state is lowered by 1.5 kcal/mol, similar to the predic- 
tions of Ref. 31. 



B. Choice of exchange correlation functional 

While the BLYP exchange correlation functional used 
in this work reproduces the water-water pair correlation 
function well, 7 it is known to overestimate correlation 
effects in hydrogen bonded systems and slightly under- 
estimate proton transfer energy barriers. For example, 
from a survey of proton transfer barriers in molecular 
systems, 69 ! 70 / 1 ! 72 B3LYP and MP2 results are found to 
be within 0.6 kcal/mol with the highest level of quan- 
tum chemistry calculations (QCISD(T) and CCSD(T)) 
in a suite of test cases, while the BLYP AE* is at most 
2.3 kcal/mol less than B3LYP predictions. This suggests 
that the accuracy of the small, IIp^ZW proton trans- 
fer barrier predicted using the BLYP functional can be 
assessed by comparing with B3LYP results in a suitable 
basis set. 

Again we consider a supermolecule geometry simi- 
lar to that in Ref. 31, using the same basis and re- 
fined basis sets as before The B3LYP functional pre- 
dicts the aqueous phase IIp^ZW AG and AG* are pre- 
dicted to be -4.0 kcal/mol and 2.6 kcal/mol, respec- 
tively. These predictions are similar to Kassab et al.'s 
AG = -5.4 kcal/mol and AG* = 2.2 kcal/mol, although 
those were obtained using a smaller basis set. 

The BLYP functional, which we use in our AIMD sim- 
ulations, predicts that the aqueous phase Hp— >ZW AG 
and AG* are -5.4 kcal/mol and 1.2 kcal/mol, respec- 



tively. This confirms that calculations with the BLYP 
functional yield similar but slightly smaller AG*'s than 
with the B3LYP functional for the proton transfer sys- 
tems of interest. The more reliable B3LYP model will 
raise AG* slightly, but this is opposite in sign to the ZPE 
correction and the two partially cancel. The BLYP and 
B3LYP AG's also differ by only 1 kcal/mol. Yet the AG 
of —5.4 kcal/mol predicted using the BLYP functional 
with a supermolecule plus dielectric continuum approach 
is considerably smaller than the AIMD umbrella sam- 
pling prediction of -11.2 kcal/mol, even after the latter is 
corrected for ZPE. 

Thus there is a several kcal/mol discrepancy between 
our BLYP-based AIMD results, and what we compute us- 
ing BLYP and a static 3- water supermolecule plus dielec- 
tric continuum calculation. We tentatively assign this to 
the fact that 3 water molecules are not sufficient to model 
the -NELj" and -COO - first hydration shells, although the 
supermolecule predictions also show some dependence on 
basis set and choice of dielectric continuum model. The 
limitations of using a small number of water molecules 
have been pointed out already in the literaturei 73 ! 74 



Appendix B: Electrostatic boundary conditions and 
finite size effects in AIMD simulations 

In this section we describe the finite size corrections to 
the hydration free energy AGhyd of a glycine zwitterion 
in water, which is closely related to the potential of mean 
force AG(R) associated with the intramolecular proton 
transfer in glycine. We consider two estimates of this 
correction: that due to a dielectric continuum argument, 
and explicit calculations of the AIMD potential of mean 
force profile by varying the simulation box size. Both 
indicate that the 11.8 A simulation box is adequate for 
glycine zwitterion in water. 

The following arguments are general and illustrate the 
ability of ab initio molecular dynamics (AIMD) meth- 
ods to successfully predict the hydration free energies of 
dipolar species in water. 



Dielectric Continuum Estimate 

The effect of periodic boundary conditions on liq- 
uid state computer simulations has been a well-studied 
subject 2^ In our case, zwitterion glycine has a relatively 
large dipole moment; we compute a value of order 15 
Debye in the gas phase. Nevertheless, using a dielec- 
tric continuum estimate, we will show that the relatively 
small box size has little effect on the glycine zwitterion 

AGhyd- 

Ewald sums and periodic boundary conditions are al- 
most invariably used in AIMD simulations, including 
those reported in this work. Consider the AGhyd of a sin- 
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size contribution -Edipoie in Eq. [21 is replaced by 
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FIG. 8: Illustrations of dielectric continuum estimates for 
a dipolar solute embedded in a high dielectric solvent. Hol- 
low triangle: solute; filled triangle: solvent molecules. The 
triangles point in the directions of the dipole moments. The 
black circles illustrate which interactions between the solute 
and the periodic images are accounted for. (a) solute-solute 
image interaction only (attractive, corrected by Eq.[3Jl; (b) all 
interactions (repulsive, corrected by Eq. |3J ; (b) only solute- 
solvent image (repulsive, corrected by Eq. 2J. With pair- wise 
additive force fields, (a) and (c) add up to (b); thus they 
cancel to a large extent in a medium with a high dielectric 
constant. Despite the particulate representation of solvent 
depicted herein, for the purpose of this discussion, the sol- 
vent is treated as a dielectric continuum. 



gle dipolar solute (molecule) at infinite dilution computed 
using Ewald summation. In a vacuum (Fig. 8a), the 
solute-solute image contribution to the computed con- 
figurational energy is given by2& 



solute— solute 



-(27rd 2 /3i 3 



(2) 



in atomic units, where d is the dipole moment and L 
is length of the cubic simulation cell. This solute-solute 
image interaction due to the Ewald sum is attractive in 
a cubic box; the boundary conditions overestimate the 
magnitude of AGhyd- The correct energy of a dipolar 
solute in vacuum is then the value computed using Ewald 

Sum minus Solute- solute- 

For a box size of L = 11. 8A and |d| = 15 Debye, Eq.[3] 
implies a ~ +5 kcal/mol correction in vacuum. 

Our zwitterion glycine molecule, however, resides in 
water, not in vacuum (Fig. 8b). In this case, the finite 



E. 



screened 



-(2^d 2 /3£ 3 )/e c 



(3) 



where e is the (relative) static dielectric constant of 
water. Therefore, the correct result is the Ewald sum 
value plus |-E scree ned|- A similar dielectric continuum es- 
timate of screened solute-solute image interaction has 
been used successfully to understand the finite size ef- 
fect in the potential of mean force of sodium chloride ion 
pair separation, computed using force field based molec- 
ular dynamics. 77 In that case, as the sodium and chloride 
ions are pulled apart from their contact ion-pair configu- 
ration, a large dipole moment is incurred. Nevertheless, 
due to the strong aqueous dielectric screening, the use of 
Ewald sum and periodic boundary conditions still lead 
to results well converged with box size (see Figs. 1 & 2 
of Ref. 77.) Estimates similar to Eq. [3] are also used in 
solid state density functional theory calculations of defect 
energetics 3& 

Assuming AIMD water exhibits e ~ 80, the dielectric 
continuum estimate of Eq.[3] suggests that the Ewald sum 
we use entails a correction of less than +0.1 kcal/mol to 
AGhyd for our box size of L = 11.8 A. 

In hydration free energy calculations that employ clas- 
sical force fields, the corrections due to Ewald sums and 
periodic boundary conditions are often described from a 
different perspective^ Unlike the case with density func- 
tional theory, where the interaction is inherently many- 
body and not pairwise decomposible, classical force fields 
typically allow the computation of strictly solute-water 
interactions. Thus the interaction of a dipolar solute with 
all periodic images, Fig. 8b, can be unambiguously sepa- 
rated into solute-solute image (Fig. 8a) and solute-water 
image (Fig. 8c) contributions. 

Using explicit force field-based molecular dynamics 
simulations, Hummer, Pratt, and Garcia*£ have elo- 
quently discussed the ramification of this separation. If 
only the solute- water terms is used, the long range solute- 
water image interaction, pictorially depicted in Fig. 8c, 
can lead to considerable finite size dependence in AGhyd- 
If the entire expression of Fig. 8b is used, i.e., if solute- 
solute self-energy term (Fig. 8a) is also included, fi- 
nite size dependences become negligableS From such an 
analysis, it can be deduced that the solute-water image 
interaction is repulsive: 



E. 



solute- 



-water « +(2^d 2 /3L 3 ), 



(4) 



Equation 0] decreases AGhyd- It is equal and opposite to 
Eq. [21 the self-energy Makov and Payne corrected for a 
dipole in vacuum, as it should be; compare the expres- 
sions in Eq. 19, Ref. 78, with the last term of Eq. 6 in 
Ref. 76. 

The physics behind the success of the solute-solute 
image self-term in removing finite size effects from the 
solute-water image interaction (Fig. 8c) is as follows. 
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The solute with dipole moment d causes water dipole 
moments in its vicinity to align against it. Since water 
has a large dielectric constant, to a zeroth approxima- 
tion the screening of the solute dipole moment is com- 
plete, and water molecules proximal to the solute exert 
a ~ d dipole moment per simulation cell. (On aver- 
age, the simulation cell should have a net dipole moment 
that approaches zero. For an anology, consider a metal, 
which has an infinite dielectric constant and cannot sup- 
port internal electric fields; in that case, it is obvious each 
unit cell has a zero dipole moment, as long as localized 
Wannier functions are judiciously used to demarcate the 
centers of charge.) Thus the solute- water images in the 
periodically replicated system give a repulsive energy of 
+ (27rd 2 /3L 3 ) (Fig. 8b), equal and opposite to the solute- 
solute image interaction. 

The solute-solute image and the solute- water image in- 
teractions are not pairwise additive in AIMD simulations, 
where energies are inherently many-body in nature. In- 
stead, AIMD propagates Newton's equations of motion 
according to the total energy and forces computed us- 
ing Ewald summation, which automatically include both 
long-range solute-solute image (Fig. 8a) and solute-water 
image interactions (Fig. 8c) . Thus AIMD simulations do 
not suffer from the large finite size correction associated 
with Eq. 21 namely the unfavorable solute-water image 
interaction (Fig. 8c) — as long as e is large and the 
dielectric continuum approximation is valid. 

Since water does not have an infinite dielectric con- 
stant, the screening of the solute dipole is incomplete, 
and the solute-water image repulsion (Fig. 8c) should 
be slightly less than the solute-solute image attraction 
(Fig. 8a) . Equation [3] is an estimate of this residual ef- 
fect. In the t — > oo limit, the system is metallic and 
this correction vanishes. When e Q = 1, the system is in a 
vacuum, and Eq. Incorrectly reduces to Eq.|3 For e « 2, 
the above arguments may not be valid because they rely 
on the assumption of strong dielectric screening. 

In summary our dielectric continuum estimate indi- 
cates that the finite size correction for our AIMD glycine 
zwitterion hydration free energy is less than 0.1 kcal/mol. 



Explicit Simulation Results 

The above analysis assumes a dielectric continuum de- 
scription of water. The particulate nature of water is not 
taken into account; with a relatively small simulation cell 
size and a small number of water molecules, it may not be 
completely valid. To investigate this, we have explicitly 
varied the simulation cell size in AIMD calculations of the 
glycine intramolecular proton transfer potential of mean 
force (AG(R)) in two umbrella sampling windows. We 
use a simulation cell roughly twice the volume used in the 
main text, containing 98 instead of 52 water molecules. 
10 ps AIMD trajectories are used for collecting statistics. 
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FIG. 9: Finite size effect in AG(R) for glycine in sampling 
window (f). The blue lines demarcate the R range used in this 
window for the overall potential of mean force. This is the 
unconstrained zwitterion window; a wide range of R values 
are sampled, and the statistics are not as good as in other 
windows. 

First we consider the unconstrained zwitterion sam- 
pling window (a). See Fig. 9. This sampling window does 
not use constraints, and so the glycine zwitterion con- 
formations fluctuate more than in other windows. As a 
result, the statistics for AG(R) tend to be worse than for 
other windows. Nevertheless, it can be seen from Fig. 9 
that AG(R) predicted for the two simulation cell sizes 
are within statistical uncertainties of each other within a 
range of R where AG(R) does not vary by more than 3 
ke T. The mean hydration structures are similar in both 
trajectories as well. The average hydration numbers of 
the larger and smaller simulation cells are with 0.1 wa- 
ter molecules of each other, with the former having 20% 
more intramolecular hydrogen bond — well within the 
statistical uncertainty. 

Since the AG(R) profile in window (a) is relatively flat, 
we also look at window (d). This window contributes 
half (~ 5 kcal/mol) of the overall AG(R), and should 
shed more light on finite size effects. In Fig. 10, we see 
that the larger and smaller simulation cells yield AAG of 
4.5 and 5 kcal/mol, respectively. The difference is within 
the statistical uncertainty. From these explicit simulation 
results, we conclude that the finite size effect should be 
within our statistical uncertainty. 
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FIG. 10: Finite size effect in AG(R) for glycine in sampling 
window (d). The blue lines demarcate the R range used in 
this window for the overall potential of mean force. 
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